Mathematical modelling of vaccination rollout and NPIs lifting on COVID-19 transmission with VOC: a case study in Toronto, Canada

Background Since December 2020, public health agencies have implemented a variety of vaccination strategies to curb the spread of SARS-CoV-2, along with pre-existing Nonpharmaceutical Interventions (NPIs). Initial strategies focused on vaccinating the elderly to prevent hospitalizations and deaths, but with vaccines becoming available to the broader population, it became important to determine the optimal strategy to enable the safe lifting of NPIs while avoiding virus resurgence. Methods We extended the classic deterministic SIR compartmental disease-transmission model to simulate the lifting of NPIs under different vaccine rollout scenarios. Using case and vaccination data from Toronto, Canada between December 28, 2020, and May 19, 2021, we estimated transmission throughout past stages of NPI escalation/relaxation to compare the impact of lifting NPIs on different dates on cases, hospitalizations, and deaths, given varying degrees of vaccine coverages by 20-year age groups, accounting for waning immunity. Results We found that, once coverage among the elderly is high enough (80% with at least one dose), the main age groups to target are 20–39 and 40–59 years, wherein first-dose coverage of at least 70% by mid-June 2021 is needed to minimize the possibility of resurgence if NPIs are to be lifted in the summer. While a resurgence was observed for every scenario of NPI lifting, we also found that under an optimistic vaccination coverage (70% coverage by mid-June, along with postponing reopening from August 2021 to September 2021) can reduce case counts and severe outcomes by roughly 57% by December 31, 2021. Conclusions Our results suggest that focusing the vaccination strategy on the working-age population can curb the spread of SARS-CoV-2. However, even with high vaccination coverage in adults, increasing contacts and easing protective personal behaviours is not advisable since a resurgence is expected to occur, especially with an earlier reopening. Supplementary Information The online version contains supplementary material available at 10.1186/s12889-022-13597-9.


Model
The system of ODEs describing the dynamic is given by: The list of variables and assumptions is given in Table SI1. Only susceptible individuals, aged 10 years and older, will receive the vaccine 2. Immunity follows two steps: partial (receiving one dose) and full (receiving two doses) 3. The vaccine efficacy is age-dependent (higher for teenagers and adults, lower for elderly) 4. The vaccine efficacy is the same against wildtype variant and VOC 5. The second dose is given after 112 days (in some predictive scenarios after 50 or 21 days), following the suggestion announced by the Government of Ontario in March 2021 6. Immunity wanes from one dose of vaccine after 120 days and from two doses after 365 days 7. We assume that the coverages in Table 2

Proportion of VOC cases
To capture the increasing trend of cases from VOC, we defined a time-dependent function ( ( )) following a sigmoid function. Fig. A1 shows the proportion of cases from VOC from data (red circles) and the function used to reproduce their trend (blue curve). According to data up to May 19, 2021 the proportion of cases from VOC in Toronto reached a maximum of 0.8 by May 11, 2021. Hence, we consider 80% to be the maximum of cases generated by the new variant.

Data fitting
To calibrate the model's parameters, we employed cumulative and daily cases and deaths, and hospitalizations ( Figure SI2). Using the Latin Hypercube Sampling (LHS), we generated 500 samples for the initial guess of each parameter using a normal distribution. Then, for each initial guess of parameter set, employed the fmincon function in the MATLAB optimization Toolbox 14 to find the local minimum of the sum of squared differences between observed data and the model's estimates of daily confirmed cases and deaths, cumulative cases and deaths and hospitalizations. After finding the best parameter set for each sample, we evaluated the mean value and the standard deviation, obtaining the confidence interval where our parameters lie.

Permutations of model's analysis
All the scenarios used for the projections are shown in Figure SI3. Each scenario is described by taking one element in each column. Figure SI3: Outward-facing model coverages and base line for model's analysis. All these coverages are reached by June 14, 2021. In brackets, we report the daily doses. Each scenario is described by taking one element in each column. Figure SI4: Variation of hospitalizations with respect to the parameters estimated in the confidence interval

Contact matrix
We used the total contact matrix from a recent Canadian study 3 . However, the age groups used in this study were defined by a 5-year band from 0 to 80+. Our model is using larger age groups, then it was necessary to aggregate the original contact matrix in less groups.
Let's define the population size of age group ∈ {1,2,3, … 17}, where 1 = 0 − 4 , 2 = 5 − 9 , … , 17 = 80 + years. To better approximate the contact rates, we calculated, from the original 17x17 matrix ( ), the total contacts that an age group has with all the other age groups. To obtain this, we multiplied all the age groups by their own population size, i.e.
× . Then, to aggregate some age groups, we averaged the total contacts as follows: • For same ages belonging to new aggregation: we summed up the diagonal entries of the submatrix related to the age groups to aggregate and the average of the mixed contacts (̂= ∑ + ∑ + 2 ). For example, the new contact of the aggregated group 0-9, given by group 1 and 2, will be 11 + 22 + 12 + 21 Once we reduced the total contacts into a smaller matrix, we re-parametrized each entry of the new age group dividing the obtained contacts by the population size of the aggregate age group (i.e., = / ∑ ). Table SI3 represents the compacted matrix.  Figure SI5: Contour plots of Rc assuming that the following coverages reached for age groups 10-19, 60-79 and 80+ years are 20%, 80%, and 90%, respectively, when the NPIs level reopening is (A) none, (B) partial, (D) total and (E) pre-pandemic. As expected, as the vaccination coverage increases, the values of the reproduction number decrease. Also, we observe that with the lowest reopening level, to reduce the reproduction number below 1, it is sufficient to vaccinate age groups 20-39 and 40-59 years above 60% and 62%, respectively. On the other hand, a relaxation of NPIs and increase in contacts as in NPIs partial reopening, the Rc will always be greater than 1. Similar results, but higher , are shown with NPIs pre-pandemic reopening (C). A B C D Table SI4: Percentage change of cumulative cases and deaths with respect to the base line NPIs no reopening in SI Figure SI3 with partial reopening in September, when age groups 60-79 and 80+ reached coverages 80%, 90% by June 14. Cases and deaths are reported comparing different coverages for age group 10-19 years, assuming 40-59 years fixed at 70% coverage (top table) Table SI5: Percentage change of cumulative cases with respect to the baseline NPIs no reopening in SI Figure  SI3 with partial, total and pre-pandemic reopening in August and September, when age groups 10-19, 60-79 and 80+ reached coverages 20%, 80%, 90%. The second dose is given at a rate of 1/112 days -1 .

Identification of the best combination of vaccination coverages and NPIs lift dates
Projected percentage change of cumulative cases with respect to the base line NPIs no reopening in SI Figure SI3 In reopen in AUGUST   Figure SI3 with partial, total and pre-pandemic reopening in August and September, when age groups 10-19, 60-79 and 80+ reached coverages 20%, 80%, 90%. The second dose is given at a rate of 1/112 days -1 .

Effect of reducing time between first and second dose
Projected percentage change of cumulative cases with respect to the base line NPIs no reopening in SI Figure SI3 with Table SI10: Percentage change of cumulative deaths with respect to the base line NPIs no reopening in SI Figure SI3 with partial, total and pre-pandemic reopening in September and second dose given after 21 or 50 days. Age groups 10-19, 60-79 and 80+ are assumed to reach coverages 20%, 80%, 90% by mid June. Par.= partial; Tot.= total; Pre-pan.= pre-pandemic.
Projected percentage change of cumulative deaths with respect to the base line NPIs no reopening in SI Figure SI3

Sensitivity Analysis
Using the Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) we conducted sensitivity analysis on the parameters related to vaccination as well as infection-related parameters.  Table SI11 shows the PRCCs of the sampled parameters , ∈ {2,3,4,5,6}, and , the daily doses in age group i, and the rate of receiving the second dose, respectively, on the cumulative cases and deaths. We observe that the age groups 3 and 4, namely, 20-39 and 40-59 years present the highest PRCC among the daily doses, suggesting that an increased vaccine coverage of these age groups leads to the largest reduction in cases and deaths. Moreover, is negatively correlated to cases and deaths, suggesting that if this rate is small, hence the time between doses is longer, cases and deaths will increase. Similar results are visible for the hospitalizations reported 50 days after reopening in June. Table SI12 shows the PRCC of some of the infection-related parameters on the model outcomes. Increase of contact, susceptibility of adults aged between 20 and 59 years show a significant positive correlation on deaths and cases, suggesting that reopening stages and higher susceptibility of adults will generate an increase of the infection, Cases data until December 2021 Figure SI9: Daily cases reported in Toronto from December 2020 to December 2021 15 (green). The dashed lined represents the end of the period of time used to calibrate the model. The curves represent the model outcomes, with highest coverage among adults, under the following scenarios: total reopening in August with second dose given 21 days after first dose; total reopening in September with second dose given 112 days after first dose and lower efficacy; total reopening in September with second dose given 21 days after first dose; partial reopening in September with second dose given 112 days after first dose.
Data following the period of time used for the model calibration show similar trend to our model prediction, with a decrease trend until August 15, followed by a slight increase. towards the end of 2021, we observe a sharp increase, attributable to the emergence of Omicron.